library(sf)
library(rgee)
ee_Initialize('denis.berroeta@gmail.com', drive = TRUE)7 S5: Índices Espectrales
Índices Espectrales con R y Rgee
7.1 Introducción
Los índices espectrales son herramientas utilizadas en percepción remota para analizar la reflectancia de la superficie terrestre en diferentes longitudes de onda del espectro electromagnético. Estos índices se basan en la relación entre la reflectancia en dos o más bandas espectrales y se utilizan para medir diferentes propiedades de la vegetación, como la salud, la humedad y la densidad.
Algunos de los índices espectrales más utilizados son:
- NDVI (Índice de Vegetación de Diferencia Normalizada): Este índice se utiliza para medir la salud y densidad de la vegetación. Se calcula utilizando la reflectancia en la banda roja y la banda infrarroja cercana. Un valor alto de NDVI indica una vegetación densa y saludable.
- EVI (Índice de Vegetación Mejorado): Este índice también se utiliza para medir la salud y densidad de la vegetación, pero tiene en cuenta factores como la influencia del suelo y la atmósfera. Se calcula utilizando la reflectancia en la banda azul, roja e infrarroja cercana.
- SAVI (Índice de Vegetación Ajustado para Suelo): Este índice se utiliza para reducir la influencia del suelo en la medición de la vegetación. Se calcula utilizando la reflectancia en la banda roja y la infrarroja cercana y tiene en cuenta la cantidad de luz reflejada por el suelo.
- NLI (Índice de Línea de Agua Normalizado): Este índice se utiliza para medir la cantidad de agua en una superficie. Se calcula utilizando la reflectancia en la banda verde y la banda infrarroja cercana.
- NDWI (Índice de Agua Normalizado): Este índice también se utiliza para medir la cantidad de agua en una superficie, pero tiene en cuenta la reflectancia en la banda del infrarrojo de onda corta y la banda verde.
Estos índices espectrales son solo algunos ejemplos de las muchas herramientas disponibles para analizar la vegetación y el agua utilizando percepción remota.
7.2 Índices Urbanos
Algunos ejemplos de índices espectrales normalizados utilizados para analizar el entorno urbano con percepción remota son:
- NDBI (Índice de Brillo de la Edificación Normalizado): Este índice se utiliza para identificar y mapear áreas urbanas construidas. Se calcula utilizando la reflectancia en la banda del infrarrojo cercano y la banda del rojo.
- SI (Índice de Suelo): Este índice se utiliza para detectar y mapear la presencia de suelo desnudo en áreas urbanas. Se calcula utilizando la reflectancia en la banda del rojo y la banda del infrarrojo cercano.
Estos índices espectrales normalizados pueden ser útiles para el análisis de la planificación urbana, la gestión del medio ambiente y otros fines relacionados con la ciudad.
Haciendo uso de proyecto awesome spectral indices que ofrece a los usuarios las herramientas para consultar índices espectrales y las bandas requeridas para el cálculo de un índice. Incluso cuentan con Aplicación web
7.3 Calculo de Índices de Urbanos
Cargar Librerías
library(tidyverse)
library(ggplot2)
library(viridis)
library(mapview)Si seleccionamos los índices espectrales obtendremos la siguente lista
Seleccionares índices espectrales urbanos (comparaciones de indices urbanos) y escogeremos el NDBI (Normalized Difference Built-Up Index).
7.3.1 Selección de Área de Estudio
mi_comuna <- "LAS CONDES"
roi <- readRDS("data/rds/zonas_urb_consolidadas.rds") %>%
st_as_sf() %>%
filter( NOM_COMUNA == mi_comuna) %>%
st_transform(4326)
plot(roi$geometry)
roi_ee <- roi%>%
st_union() %>%
st_geometry() %>%
sf_as_ee()7.3.2 Selección Imagen Satelital
Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A
maskS2clouds <- function(image) {
qa <- image$select("QA60")
# Bits 10 and 11 are clouds and cirrus, respectively.
cloudBitMask <- bitwShiftL(1, 10)
cirrusBitMask <- bitwShiftL(1, 11)
# Both flags should be set to zero, indicating clear conditions.
mask <- qa$bitwiseAnd(cloudBitMask)$eq(0)$And(
qa$bitwiseAnd(cirrusBitMask)$eq(0)
)
# Return the masked and scaled data, without the QA bands.
image$updateMask(mask)$
divide(10000)$
select("B.*")$
copyProperties(image, list("system:time_start"))
}# Map the function over one year of data and take the median.
# Load Sentinel-2 TOA reflectance data.
S2 <- ee$ImageCollection("COPERNICUS/S2_HARMONIZED")$
filterDate("2016-01-01", "2016-12-31")$
filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", 20))$
map(maskS2clouds)$
median()$
clip(roi_ee)viz <- list(bands = c("B4", "B3", "B2"), min = 0, max = 0.3)
Map$centerObject(eeObject = S2,zoom = 12)
Map$addLayer(S2, viz, "RGB")7.3.3 Aplicar Índice
la función NDBI corresponde a la siguiente:
NDBI = \frac{S1-N}{S1+N} Por lo anterior es muy importante conocer los nombres de la bandas para crear correctamente la función, a continuación de exponen las bandas del producto satelital sentinel.
F_NDBI<- function (image){
ndbi<-image$expression('float(S1-N)/float(S1+N)',
opt_map=list("N"=image$select("B8"),
"S1"=image$select("B11")))$
rename("NDBI")
return(ndbi)
}S2_NDBI <- F_NDBI(S2)7.3.4 Visualizar Resultados
pal_ndbi <- c('white' , 'white', 'white', 'red' , 'red', 'red')
viz_ndbi<- list(min=-1,
max=1,
palette=pal_ndbi)
Map$centerObject(eeObject = roi_ee,zoom = 12)
Map$addLayer(S2_NDBI,visParams =viz_ndbi)7.3.5 Calcular un índice
mz_sii <- readRDS("data/rds/mz_constru_SII.rds") %>%
st_as_sf() %>%
filter(n_com == mi_comuna) %>%
mutate(id = 1:nrow(.))Extraer los valores con la función ee_extract cuya documentación se encuentra en el siguiente link
#Extract values - getInfo
ee_com_urb <- ee_extract(
x = S2_NDBI,
y = mz_sii["id"],
scale = 10,
fun = ee$Reducer$mean(),
sf = TRUE
)mz_sii_urb <- mz_sii %>%
left_join(ee_com_urb %>% st_drop_geometry(), by = "id")7.3.6 Corrección
mz_sii_urb <- mz_sii_urb %>%
mutate(NDBI = ifelse(total == 0, NA, NDBI))7.3.7 Histograma de NDVI
ggplot(mz_sii_urb, aes(x = NDBI)) +
geom_histogram(bins = 50, fill = "#1d91c0", color ="gray80")+
theme_bw()+
labs(title="Histograma de Índice de Urbanismo", x ="NDBI", y = "Frecuencias")+
theme(plot.title = element_text(face = "bold",colour= "gray60", size=10)) 
7.3.8 Visualización espacial de NDBI
mapview(mz_sii_urb, zcol= "NDBI")library(ggplot2)
library(viridis)
# Visualización ggplot y sf
ggplot() +
geom_sf(data = mz_sii_urb, aes(fill = NDBI), alpha=0.7)+
scale_fill_viridis_c(option = "E", direction = -1)+
ggtitle("Urban Index: NDBI") +
theme_bw() +
theme(legend.position="none")+
theme(panel.grid.major = element_line(colour = "gray80"),
panel.grid.minor = element_line(colour = "gray80"))
7.3.9 Guardar los resultados
saveRDS(ee_com_urb, file = "data/rds/ndbi_2016.rds")7.4 Monitoreo de Sequía Laguna de Aculeo
Para el siguiente análisis utilizaremos el NDWI: Normalized Difference Water Index)
7.4.1 El Problema
Noticias Relacionadas:
“La sequía de Aculeo vista desde el espacio”. Diario La Tercera 11-MAY-2018
7.4.2 Lectura de Insumo Cuerpos de Agua Nacional
masas_agua <- st_read("data/shape/masas_agua.shp", quiet = T)
masas_aguaSimple feature collection with 16644 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -5659.382 ymin: 3803587 xmax: 692350 ymax: 8050958
Projected CRS: WGS 84 / UTM zone 19S
First 10 features:
NOMBRE TIPO_MAGUA COMUNA AREA_KM2 RULEID
1 Laguna Tebinquiche Laguna <NA> 2.08693247 5
2 <NA> Laguna <NA> 0.84982790 5
3 Salar de Pujsa Salar <NA> 0.22226662 7
4 Salar de Pujsa Salar <NA> 16.64594244 7
5 Salar de Atacama Salar <NA> 49.61221690 7
6 Salar de Pujsa Salar <NA> 0.02636444 7
7 Salar de Loyoques o Quisquiro Salar <NA> 65.49493327 7
8 Salar de Loyoques o Quisquiro Salar <NA> 0.63644089 7
9 Salar de Pujsa Salar <NA> 0.08644732 7
10 Laguna en Salar de Loyoques o Quisquiro Laguna <NA> 0.02395804 5
SHAPE_Leng SHAPE_Area geometry
1 12353.4677 2086932.47 MULTIPOLYGON (((575205.1 74...
2 4397.5639 849827.90 MULTIPOLYGON (((661560.1 74...
3 3305.8444 222266.62 MULTIPOLYGON (((652931.2 74...
4 26208.2874 16645942.44 MULTIPOLYGON (((652560 7434...
5 44044.4984 49612216.90 MULTIPOLYGON (((569340.4 74...
6 606.6370 26364.45 MULTIPOLYGON (((653547.4 74...
7 67250.4947 65494933.27 MULTIPOLYGON (((672122.6 74...
8 5459.5505 636440.89 MULTIPOLYGON (((671993 7433...
9 1147.0949 86447.32 MULTIPOLYGON (((653825.7 74...
10 686.3645 23958.04 MULTIPOLYGON (((673224.5 74...#Buscar nombre del cuerpo de aguar por patrón de texto
grep(pattern = "Aculeo$", x = masas_agua$NOMBRE, value = T)[1] "Laguna de Aculeo"# Leer Laguna Aculeo
lag_aculeo <- masas_agua%>%
filter(NOMBRE == "Laguna de Aculeo")%>%
sfheaders::sf_remove_holes() ## remover poligonos internos%>%>%plot(lag_aculeo$geometry,col = "skyblue", main = "Laguna Aculeo")
#Laguna Formato ee
lag_aculeo_ee <- lag_aculeo%>%
st_buffer(500)%>%
st_transform(4326)%>% #crs atlon
sf_as_ee()
# regiòn para definir regiòn de interès para los mapas
region <- lag_aculeo_ee$geometry()$bounds()7.4.3 Lectura de insumo satelital visualización previa
imagen <- ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
filterBounds(region)$
filterDate('2018-01-01','2018-12-31')$
filterMetadata('CLOUD_COVER','less_than', 5)$
median()$
clip(region)
ndwi <- imagen$normalizedDifference(c("B3", "B5"))$clip(region)
ndwiMask <- ndwi$updateMask(ndwi$gte(0))#greater than or equal to the given value.
ndwiViz <- list(min = 0, max = 1, palette = c("white", "#00FFFF",'#0080FF', "#0000FF"))
imgViz <- list(min = 0, max = 0.5, bands = c("B4", "B3", "B2"), gamma = c(0.95, 1.1, 1))
Map$centerObject(region, zoom = 13)
mapa <- Map$addLayer(imagen, imgViz,'RGB')+
Map$addLayer(ndwiMask, ndwiViz,'NDWI')
mapa7.4.4 Construcción de un Mosaico con NDWI y Falso Color
imageRGB <- imagen$visualize(bands = c('B5', 'B4', 'B3'), max = 0.5, gamma = c(0.95, 1.1, 1))
ndwiRGB <- ndwiMask$visualize(min = 0, max = 1, palette = c('#00FFFF', '#0000FF'))
mosaic <- ee$ImageCollection(c(imageRGB, ndwiRGB))$mosaic() #función crea mo
Map$centerObject(region, zoom = 13)
Map$addLayer(mosaic , list(), 'mosaic')7.4.5 Creación de Serie Anual desde el 2015 al 2020
imagen <- ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
filterBounds(region)$
filterDate('2015-01-01','2020-12-31')$
filterMetadata('CLOUD_COVER','less_than', 5)
anual <- ee$List$sequence(2015, 2020)
anual_l8 <- function(y) {
imagen$filter(ee$Filter$calendarRange(y, y, "year"))$
median()$
select("B4", "B3", "B2")$
clip(region)}
l8_year <- anual$map(ee_utils_pyfunc(anual_l8))
l8_year_2015 <- ee$Image(l8_year$get(0))
l8_year_2016 <- ee$Image(l8_year$get(1))
l8_year_2017 <- ee$Image(l8_year$get(2))
l8_year_2018 <- ee$Image(l8_year$get(3))
l8_year_2019 <- ee$Image(l8_year$get(4))
l8_year_2020 <- ee$Image(l8_year$get(5))
visparams <- list(
bands = c("B4", "B3", "B2"), min = 0,max = 0.6,gamma = 1.4)
Map$centerObject(region, zoom = 13)
map1 <- Map$addLayer(l8_year_2015, visparams, name = "2015")+
Map$addLayer(l8_year_2016, visparams, name = "2016")+
Map$addLayer(l8_year_2017, visparams, name = "2017")+
Map$addLayer(l8_year_2018, visparams, name = "2018")+
Map$addLayer(l8_year_2019, visparams, name = "2019")+
Map$addLayer(l8_year_2019, visparams, name = "2020")
map1# referencias: # https://github.com/r-spatial/rgee/blob/examples/ImageCollection/creating_monthly_imagery.R